function model = createToyModelForgapFind()
% Create a toy model that has a partially disconnected cycle (R4, R5, R6), a missing exporter(For I) and
% a missing importer (For H)


model = createModel();
%Reactions in {Rxn Name, MetaboliteNames, Stoichiometric coefficient} format
Reactions = {'R1',{'A','B'},[-1 1],0,1000;...
             'R2',{'B','C'},[-1 2],0,1000;...
             'R3',{'C','D'},[-1 1],0,1000;...
             'R4',{'E','F'},[-1 1],0,1000;...
             'R5',{'F','G'},[-1 1],0,1000;...
             'R6',{'G','E'},[-1 1],0,1000;...
             'R7',{'G','D'},[-1 1],0,1000;...
             'R8',{'H','C'},[-1 1],0,1000;...
             'R9',{'C','I'},[-1 1],0,1000;...
             'EX_A',{'A'},[-1],-1000,1000;...
             'EX_D',{'D'},[-1],-1000,1000};
ExchangedMets = {'A','D'};
%Add Reactions
for i = 1:size(Reactions,1)
    %All reactions are irreversible
    model = addReaction(model,Reactions{i,1},'metaboliteList',Reactions{i,2},...
          'stoichCoeffList',Reactions{i,3},'lowerBound',Reactions{i,4},'upperBound',Reactions{i,5},...
          'printLevel',-1);
end

%Add Exchangers
%model = addExchangeRxn(model,ExchangedMets,-1000*ones(numel(ExchangedMets),1),1000*ones(numel(ExchangedMets),1));
model = changeObjective(model,'EX_D',1);